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Abstract 

We study the metal-insulator transition and magnetic ordering in the Hub- 
bard model using the particle-hole mapping. The analysis simplifies near the 
ferromagnetic limit. We find that the two dimensional(2D) Hubbard model 
has a charge excitation gap at half-filling for any finite U in this region on 
both the bipartite square lattice and the nonbipartite triangular lattice. In 
some cases, the system goes through a first-order phase transition to become 
a paramagnetic metal as S z is lowered. We also discuss the extension to the 
doped case. We find that in the large U limit, a single doped hole has a 

bandwidth of order of J rather than t at S z = 0. 
PACS numbers: 71.30.+h, 75.10.Jm, 05.30.Fk 
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I. INTRODUCTION 



The Hubbard model is one of the most extensively studied strongly correlated Fermion 
models in condensed matter physics. Despite its simple appearance, it has been used 
to explain a wide variety of physical phenomena such as ferromagnetism, || antiferromag- 
netism, M metal-insulator transition, j[JR| and more recently, it has been proposed as a 
model for high-T c superconductivity One of the interesting features of the Hubbard 
model is the possibility of transforming from the repulsive (U > 0) to the attractive (U < 0) 
model. It has been known for a long time that a positive-U Hubbard model can be mapped 
into a negative-U model by making a particle-hole transformation on one of the spin species. 
Here we use this mapping to study the magnetic and transport properties of Hubbard 



model on various lattices. One of the motivations to use this transformation is to make use 
of the existing knowledge about the attractive Fermion gas in finding a good trial wave- 



function for the Hubbard model. For example, Leggett [] 1] and Nozieres and Schmitt-Rink 



12] have shown that the BCS wavefunction contains the right physics of attractive Fermion 



gas in both the weak and the strong coupling limits, and it may be used as an interpola- 
tion scheme in between to describe the progressive buildup of pairing correlations in the 
ground state. Thus we can use a unified method for all U. Furthermore, this wavefunction 
(in the transformed variables) provides a natural description of the binding and unbind- 
ing of particle-hole-pairs. This binding can be regarded, in certain situations, as the Mott 
transition into the insulating state. 

Because of the importance of Hubbard model, exactly soluble limits or modifications are 
of great interest, as evidenced by much work on the Nagaoka limit of small doping and large 
U, and on the one-dimensional t-J model. By mapping the positive-U Hubbard model into 
a negative-U model, we are able to solve the Hubbard model in the limit of high total z 
spin, S z for all U. We find that dimensionality plays an important role in the metal-insulator 
transition in this limit. InlD and 2D, the Hubbard model always has a charge excitation 
gap at half-filling at high S z , and therefore, is always an insulator for any finite U in the 
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high-S^ region. In three or higher dimensions, however, there exist a critical interaction [/&, 
which separates a metallic phase (small U) from an insulating phase (large U). Using the 
information from the exact solution, we construct a BCS type of many-body wavefunction 
for the lower S z case and calculate the phase diagram in the S z -U plane. In our work, we 
choose to fix S z rather than apply an external magnetic field which acts on the spins. We 
discuss this choice in more detail below. 

Using a particle-hole transformation for spin-up electrons, the Hubbard Hamiltonian 

H=-^2 tijcl a c ja + h.c. + cl^c^cn, (1) 

ijcr i 

can be written as 

H = Ujhthj - tijdjdj + h.c. + Uj^ d \di ~Uj^ d\dih\hi. (2) 

ij ij i i 

Here h\ = and d\ = cfi are the creation operators of the "hole" and the "doublon". c\ a 
is the original electron creation operator. The "vacuum" state |0 > has an up spin electron 
at every site. In this hole-doublon representation (referred as HDR hereafter), holes and 
doublons have an on-site attractive interaction —U, and do ub Ions have a site energy U. The 
zero-radius bound state of a doublon and a hole at site i is nothing but a spin down state 



at i. In the dilute limit, the bound state pair behaves like a hard-core boson. JTj| Also note 
the change in sign of the hole kinetic energy, a crucial part of the transformation. 

The rest of paper is organized as follows: In Sec. II, we present the exact solution in the 
high total spin limit. We show, by calculating the Kohn charge stiffness constant, that the 
system is an insulator when and only when there exists a bound state solution in HDR. We 
devote Sec. Ill to the lower total spin case. In Sec. IV, we present a possible generalization of 
our many-body wavefunction at the half-filling to the hole-doped case and use it to calculate 
the minimum energy and the bandwidth of the doped hole. We summarize our results in 
Sec. V. 
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II. THE HIGH TOTAL SPIN LIMIT 



Let us start with the exactly soluble problem of a single pair of doublon and hole, i.e., 
the band is half-filled and there is only one down spin, S z = (N — 2)/2. Here N is the 
number of electrons. One of the motivation to study the single pair problem is that in 2D, 
the existence of the bound state solution is the necessary and sufficient condition of the BCS 
many-body instability, 0] which in our case corresponds to a spin excitation gap at lower 



S z , although the actual situation is more complicated because of the particle- hole mapping, 
as we shall discuss later in Sec. III. 

The system consists of a single hole and a single doublon, all other sites being occupied 
by an up spin. In this background, both particles move freely, except for their attraction 
for one another, which is zero range. This two particle problem differs from the text book 
bound state problem only in that the particles move in a band of finite width and that the 
hole prefers to sit at the top, not the bottom, of this band. The latter circumstance implies 
that the lowest energy state of a bound pair may not have zero total momentum. 

The wavefunction of the system consisting of one hole and one doublon with total mo- 
mentum K in the HDR can be written as, 

^ = E^4 /i Ld 0> - ( 3 ) 

k 

Substituting Eqs.(2) and (3) into the Schrodinger equation yields, 

AU I N U\ 

a k = ~p TTT • ( 4 ) 

h K ~ e a + € K-q 

Here e& = — J2j tifi %k rii is the dispersion relation for the non-interaction case and A is a 
normalization constant which cancels out in the eigenvalue equation. The eigenenergy of 
the system is given by 

Ir ! = -I. (5) 

We first consider the most simple case that t^ is only nonzero for nearest neighbor 
hopping on a d-dimensional hyper-cubic lattice, then = — It Y,m=i cos k m - (We assume 



the lattice constant a = 1 throughout this paper.) The doublon kinetic energy is minimum 
at k — 0, and the hole kinetic energy at k — (ir, 7T, 7r). For every given K, one can 
calculate Eg from Eq.(5). We find that the minimum value of Ek is reached at K\ = Ki = 
. . . = Kd = 7r, consistent with the kinetic energies of the particles. So Eq. (5) can be written 

as 

f{E) = I (2nY E + At E d m=l cosq m = ~ U (6) 
in the limit A" — > 00. Here E is the minimum value of Eg and the integration is over the 



first Brillouin zone. The structure of the solutions of Eq. (6) is well known. [[L4J] There is 
an extended state continuum of eigenvalue E between energy — 4td and 4td, and there may 
or may not be a bound state below the continuum lower edge E_ = —4td depending on the 
interaction U and the dimensionality d. For d=l, the integral in Eq. (6) is divergent as 1/e 
when E approaches E_ from below so that f(E — > E_) — > —00. Also f(E — > —00) = 0. 
This means that there is always a bound state solution no matter how small the interaction 
U is. d — 2 is the marginal case, the integral in Eq. (6) is logarithmically divergent as E 
approaches E_ from below. So there is a bound state for any U, although the gap between 
the bound state and the extended state continuum decreases exponentially as e~ 8lTt ^ u when 
U is small. For d=3, the integral of Eq. (6) is convergent as E approaches E_ = — 12t and 
there exists a critical such that when U < Ub, there is no bound state. is given by 

L = f*±. I (7) 

U b J (2tt) 3 4t(3-E^=i cos q m y 1 ; 

Numerical evaluation of the integral gives [/& ~ 7.9 161 Ub increases monotonically as d is 
increased past 3. 

Next we consider the single-pair hole-doublon problem on a 2D triangular lattice with 
the nearest neighbor hopping only. The sign difference of the hopping terms of doublons and 
holes due to the particle-hole transformation now becomes important. When the system is on 
a bipartite lattice, the minus sign can be simply "gauged" away by redefining h\ = s(i)c^. 
Here s(i) equals 1 or -1 when the site i belongs to two different sublattices. After this 



redefinition, the new "hole" has the same hopping term as doublons and is able to hop 
coherently with the doublon which it forms a bound state with. Now the minimum of 
energy is reached when the transformed total momentum K = 0. When the Hubbard model 
is on a nonbipartite triangular lattice, however, the sign difference of the hopping terms of 
holes and doublons in Eq. (2) can not be gauged away, i.e., the hopping terms with +t and 
—t are not symmetric any more. For example, the minimum value of the doublon hopping 
term is — 6t on a 2D triangular lattice, while the minimum value of the hole hopping term 
is only — St. So the lower edge of the extended state continuum is at E_ = —9t instead of 
— 2zt = —12t. (Here z=6 is the number of the nearest neighbors.) The eigenenergy of a 
single pair of hole and doublon on a triangular lattice is still given by Eq. (5). Now 

= — 2t[coski + cosk 2 + cos{k\ — k 2 )]. (8) 

Here we have chosen the two primitive axes of the triangular lattice as the axes of our 

— * 

coordinate system. k\ and k 2 are the components of k in the non-orthogonal basis which 
generates the reciprocal lattice. The minimum of Eg is reached when the total momentum 
of the hole-doublon pair K is at one of the corners of the hexagonal Brillouin zone, for 

— * 

example, K = (2n/3, — 2tt/3) in the non-orthogonal basis. Then Eq. (5) becomes 

A f-^L I = _I ( 9 ) 

C J {2n) 2 E + 2V3t[sin(l - qi ) + sm(f + g 2 ) + sin(f - Ql + q 2 )\ U' KJ 

Here A c = V3/2 is the area of the primitive cell and the integral is over the first Brillouin 
zone. The integral of Eq. (9) is logarithmically divergent as E approaches the lower edge 
of the extended state continuum E_ = —9t from below, which means that, just as on the 
square lattice, a single pair of hole and doublon on the nonbipartite triangular lattice also 
always forms a bound state no matter how small the interaction U is. Similarly, we can show 
that the conclusion is also valid for the 2D Hubbard model with the nearest and the next 
nearest neighbor hopping. In this case, the noninteracting energy dispersion on a square 
lattice is given by 

efc = — 2t[cosk x + cosk x ] — \t 2 cosk x cosk y . (10) 



Here the first and the second term correspond to the nearest neighbor and the next nearest 
neighbor hopping, respectively. For t > 2t 2 > 0, the minimum of is at K — (it, it), which 
is the same as the t 2 = case. We can see that the question of whether a pair of hole and 
doublon always form a bound state for any finite U is only determined by the dimensionality 
of the system. For d > 2, there exists a critical Ub, when and only when U > £4, doublons 
and holes form charge neutral bound states; While for d < 2, they always form bound states 
for any finite U. This conclusion is easily proved by considering the behavior of the integral 
in Eq.(5) near the point in k-space where — e^_£ takes on its minimum value E— This is 
a variant of standard arguments about bound states. Its interest here is that two dimension 
is more similar to one dimension than it is to three dimensions in this limit of the Hubbard 
model, namely for high total spin but for all U. This is contrary to the zero total spin case, 
where it is usually felt that one dimension is very special as regards the metal-insulator 
transition. 

The conclusion that a bound state exists for all U in one and two dimensions and above a 
critical U in three dimension is independent of lattice structures. This is well known in single 
particle quantum mechanics. We have simply reinterpreted the result as a metal-insulator 
transition at high S z . The critical dimension for this transition is two in this model. 



To verify that the bound state is indeed an insulating state, we follow Kohn [|15|] and 
calculate the charge stiffness constant D c which measures the response of the system to an 
electromagnetic field. We assume periodic boundary conditions along the x-direction and 
thread the system with a flux $, which we represent by a vector potential A = (§/N x )e x , 
where N x is the number of sites in the x-direction. The charge stiffness constant D c = 
N x d 2 E(§)/d§ 2 is a measure of how good a conductor the system is, and vanishes for an 
insulator. [p~5| , p~6| , |2Q[1 Here E(Q), the ground state energy in the presence of the flux $, is 
given by, 

Here e^($) is the single particle kinetic energy in the presence of the flux $. For a 2D square 



lattice with the nearest neighbor hopping, for example, e^($) = — 2t[cos(q x +$ /N x )+cos(q y )]. 
One can easily write down similar expressions for other cases. We expand Eq. (|11]) up to 
the second order of The nth order equation gives d n E/d& n at $ = 0. After some algebra, 
one finds that dE/d§ is always zero, and d 2 E/d& 2 is zero when and only when E < E-, 
i.e., D c = when and only when there exists a bound state below the extended state band. 



This is because that only when E < E- = —8t, one can expand Eq. (jTTj) in terms of 
<&/[E + 4t(cosq x + cosq y )}; When E > E-, this expansion is not valid in some regions of 
the momentum space and D c is finite. This conclusion also applies to the 3D cubic lattice 
and other cases. We have therefore established the one-to-one correspondence between the 
existence of the bound state and the insulating behavior of the system. The bound state 
wavefunction contains precisely the hole-doublon phase coherence necessary for insulating 
behavior. This coherence can not be achieved when hole-doublon binding is accomplished 



by Gutzwiller-type projection methods. [16 



Next we consider an anisotropic three-dimensional system with small interlayer hopping. 
This is important in the study of high-T c superconductors. As we discussed above, the 
dimensionality of the system plays a very important role in determining whether the system 
is a Mott insulator or a metal. We can expect that even a small interlayer hopping could be 
quite important because it introduces the third dimension. Mathematically, the interlayer 
hopping term provides a cutoff to the 2D logarithmic divergence of Eq. (5). In this case, 
Eq. (5) can be written as 

(1±_ I = _I (12) 

J (2tt) 3 E + At(cosq x + cosq y ) + At' cosq z U 

Here 4t'cosq z is the interlayer hopping term. As before for t' ^ 0, there exists a Ub{t r ) such 
that when and only when U > Ub(t'), Eq. (11) has a bound state solution. In Fig. 1, we 
show the critical interaction Ub(t') as a function of t'. We can see that [/& does depend on t' 
quite strongly in the small t' region. For example, at t' = O.OOlt, we still have Uj, = 2.434£ 
even though [/& = at t' = 0. The crossover to 3D behavior occurs very rapidly as t' is 
turned on. 
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When there are few doublons and holes, i.e., in the dilute limit, the overlap or interaction 
between the hole-doublon bound states (HDBS) is small and negligible, and the single pair 
HDBS solution gives a correct physical picture of the system. In ID and 2D, the Hubbard 
model always has a charge gap at half-filling at high S z , and therefore, is always an insulator 
for any finite U in the high-S" 2 regime. (At least at S z = N/2, which is obvious, and at 
S z = N/2 — 1, which we have just proved rigorously). In 3D, when the interaction U 
is smaller than Ub, there is no bound state and the system is a paramagnetic metal (or 
semimetal since the carrier density is small); when U is larger than Ub, the doublons and 
holes form charge neutral bound states and the ground state is a gas of bound doublon-hole 
pairs, a very appealing model for an insulator. [fl,|15,16 



As the number of HDBS increases, the bound states merge into an energy band which 
corresponds to the lower Hubbard band in the original electron representation. Similarly, 
the extended state continuum corresponds to the upper Hubbard band, and the gap between 
the bound state and the extended state is nothing but the Mott-Hubbard gap. In 2D, we 
have shown that a single pair of hole and doublon always form a bound state and has a gap 
no matter how small the interaction U is. The question is whether this gap will survive the 
overlap and the interaction between the HDBS pairs. That is the question we shall address 
in the next section. 

III. GENERAL TOTAL SPIN 

In this section we consider the Hubbard model at half-filling with general S z = (N — 
2M)/2, i.e., there are M spin-down electrons (or M doublon-hole pairs in HDR) in the half- 
filling N-site system. We exploit the fact that the relationship between the high spin limit 
and the lower spin case is very similar to that between the Cooper pairing problem and the 
BCS many-body problem. It is therefore natural for us to choose a BCS type of many-body 



wavefunction for the many-pair state. [14.17 
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\K, S z >= Y[(u k + v k d\h} K _ k )\U > (13) 

k 

where u\ = \{l + ^ k /E k ), v\ = ±(l-f fc /J5 fc ), El = 4A 2 + «£f, and & = e k -e K - k ~ fi d - fi h is 
the kinetic energy of the Cooper pair measured from the Fermi surface, fi d and \Xh are the 
chemical potential of the doublons and the holes, respectively, and A = S J2 k u kV k is the 
BCS gap function which is momentum-independent in this case because of the short-range 
interaction in Hubbard model. K is the total momentum of the hole-doublon "Cooper pairs" 
which can also be thought as the wavevector of the spin-density wave in the original electronic 
picture. PJTT|| This kind of pairing state was first studied by Fulde and Ferrell, and Larkin 
and Ovchinnikov in the context of superconductivity. p!7|JT8[] Obviously this wavefunction 
is not an eigenfunction of operator S z , and the labeling S z on the left-hand side of Eq.(p~3"D 
should be understood in the sense of the average expectation value of operator S z just as the 
number of particles in the BCS case. Usually K is at one of the corners of the first Brillouin 
zone as we can see from the single-pair solution. It could move away from the corners in 
some cases, which corresponds to the incommensurate spin-density wave. 

We evaluate the expectation value of the Hamiltonian Eq. (2) using the wavefunction of 



(|13|). The gap equation is then given by minimizing the energy expectation value [19 



^ V 1 = l_ , V 

N^[(e k -e K ^ k -2^ + AA^ U' 1 j 

Here 2fi = fi d + fi h . The total z spin is given by S z = N/2 — J2k v k i which, together with the 
gap equation, can be used to calculate the chemical potential \i and the gap A, 
J_V^ri e fc - e K _ k - 2/i = _ 2S^ 

NY [(6 fe -6 A -^-2/i)2 + 4A2]i/2 J N - 1 ) 

Here N is the number of sites. It is easy to demonstrate that in 2D, Eq.(^) always has a 
nonzero solution A for any finite U, which means that there is always a charge excitation 
gap in this particular wavefunction for any finite value of U. To show that, we take the limit 
N — > oo and replace the sum Y^, k with the integral over the first Brillouin zone, 

r d 2 k 1 _ 1 

C J (2vr) 2 [( efc - ei ,_,-2/i) 2 + 4A 2 ]V2 ~ Jj- ^ > 
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Here A c is the area of the primitive unit cell. A c = 1 for the square lattice and A c = \/3/2 
for the triangular lattice. The integral in Eq. (|IED is logarithmically divergent as A goes 
to zero. This is due to the fact that generally in d = 2, (er — £ K Zk ~ ^//) vanishes along a 
curve in the momentum space. That is sufficient to cause the logarithmic divergence of the 
left-hand-side of Eq. 

This reasoning would suggest that the 2D Hubbard model is an insulator at half filling 
for any finite U, as is true for the ID model, and the square lattice with the nearest neighbor 
hopping only. To understand how the reasoning breaks down, we consider three cases: 

(A) , bipartite square lattice with the nearest neighbor (NN) hopping only, 

(B) . square lattice with NN and the next nearest neighbor (NNN) hopping, and 

(C) . nonbipartite triangular lattice with nearest neighbor hopping only. 
A. Square lattice with NN hopping only 

In this case, e^_£ = — e^, and the many-body wavefunction at A = is identical 
to the noninteracting wavefunction. The fact that the gap equation always has a nonzero 
solution for any finite U means that the Fermi liquid state is unstable against the BCS state 
(|13|) at finite U. The staggered magnetization of Hubbard model is proportional to A/U. So 
obviously a nonzero A means that the system is antiferromagnetically ordered. We therefore 
conclude that the system is always an antiferromagnetic insulator for any finite U. This result 
is well known, but the present method gives a new way of visualizing it. Consider the first 
Brillouin zone in Fig. 2 (a). The doublons occupy the region near K = 0, and the holes 
occupy the region near K = (7r,7r). The energy curves are exact translates of one another 
through (7r, 7r) for all fillings, so each hole pairs with one doublon. This particular model is 
known to have perfect nesting at half filling, but note that this property is even stronger. 
There is a perfect particle-hole nesting at all values of the magnetization. As a result, the 
pairing state is always energetically favorable. 

In the large U limit, one can solve the gap equation (|TJP by expanding the gap function 
A in terms of t/U. At S z = 0, we find 



as A ->• 0. 
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and the expectation value of the Hamiltonian in the pairing state is 

E =< K, S z = 0\H\K, S z = 0>= -ANt[^j + O(^)), (18) 

which is the expected result — NJ for a Neel state, and shows that our mean-field solution 
has the correct large-U limit. Here J = 4t 2 /U. The pairing state (|T3| ) may do better when 
U is small. The reason is that in this case the size of bound states becomes large, there 
are many particles inside the region covered by a single bound state wavefunction, which 
will suppress the relative fluctuation in the interaction effects, and the BCS mean-field 
wavefunction becomes a better approximation to the exact ground state. 
B. Square lattice with NN and NNN hopping 

The noninteracting energy dispersion in this case is given by Eq. (10). The crucial 
difference between the case (A) and (B) is that for case (B), ejc-k 7^ ~ e k, and the BCS many- 
body wavefunction ( |13|) does not reduce to the noninteracting wavefunction when A = 0, 
i.e., the fact that the gap equation ([14]) always has a nonzero solution for any finite U just 
means that the wavefunction (|13|) at A = is always unstable against the finite-A BCS 
pairing state, but the A = wavefunction is not necessarily the Fermi liquid wavefunction. 
Thus in HDR, when the holes and doublons take advantage of the attractive interaction 
between them and form the BCS condensate, they pay the price of having a little higher 
kinetic energy in case (B). This can be seen in Fig. 2(b), where the doublon Fermi surface 
near K = (0, 0) can not be mapped onto the hole Fermi surface near K = (tt, it) through a 
linear transformation k' = K — k. It is easy to show that the Fermi liquid state has lower 
kinetic energy than the pairing state fll3f) . This problem becomes severe at low S z . So when 
the gain from the condensation energy is not enough to compensate the extra kinetic energy, 
the BCS state fll3|) becomes unstable to the Fermi liquid state, the system goes through a 
first order phase transition and becomes a paramagnetic metal (Fermi liquid). Once we 
have established the correct physical picture, it becomes straightforward to calculate the 
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metal-insulator phase boundary on which the energy of the many-body wavefunction ([13]) 
is equal to the energy of a true Fermi liquid state. In Fig. 3, we show the calculated phase 
diagram of the case (B) with t% = 0.25t in the S Z -XJ plane. For S z = 0, the critical interaction 
U c is about 2.3t, which is in agreement with the result of the numerical and Hartree-Fock 
calculation by Lin and Hirsch. || As ^ goes to zero, case (B) reduces to case (A), and the 
critical interaction of the metal-insulator transition U c goes to zero. 
C. Triangular lattice with NN hopping 

Case (C) is more complicated than case (B), due to the conflict between the "pairing" 
property of the wavefunction (|T3| ) and the "tripartite" character of the lattice. As a result, 
the holes of the wavefunction ([13]) only occupy three of the six corners of the hexagonal first 
Brillouin zone (BZ), leaving the low energy states of the other three corners unoccupied or 
partially occupied. The occupied three corners, as well as the unoccupied three corners, 
are related by the reciprocal lattice vectors and form two subsets, respectively, but these 
two subsets are not related by the reciprocal lattice vectors. As shown in Fig.4(a), the 
hexagonal BZ (dotted line) is equivalent to the parallelogram (solid line) which is easier to 
deal with analytically. So by choosing K to be at one of the BZ corners (or equivalently, 
two points A and B on the diagonal of the parallelogram which divided the diagonal into 
three equal segments), we force the system breaking the point group symmetry as well as 
the time reversal symmetry. In Fig. 4(a), we also show the Fermi surfaces of the doublons 
(solid line) and holes (dashed line) at S z = OAN and 0.3N, respectively. Apparently there 
are two minimum energy positions (which correspond to two distinguish sets of corners of 
the hexagonal BZ) for holes to occupy, but in the pairing state (|13|) they can only occupy 
one of them and leave the other unoccupied. This, as well as the distortion of the Fermi 
surface, costs extra kinetic energy, which opens doors for those many-body states other 
than the Fermi-liquid or the state of ( |i3|) with K at the one of the corners of BZ as the 
possible candidate of the ground state in some parameter region. One possible candidate 
is the incommensurate spin-density wave state which corresponds to a Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) type of state in HDR. In comparison, the four corners of the BZ of the 
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square lattice are related by the reciprocal lattice vector and the increasing of the kinetic 
energy in the pairing state in the NNN hoping case is solely due to the distortion of the 
Fermi surface. For this reason, we expect the critical interaction U c of the instability of 
the BCS state (|1|) in case (C) to be larger than that in case (B), and the system may go 
through other kind of states before becoming a paramagnetic metal (Fermi liquid). 

We calculate the energy of the pairing state (|13|) for general K. We find that the energy 
of the pairing states K at the corner of first Brillouin zone is at least a local minimum. It 
is quite possibly also a global minimum. Since there are two set of the optimum K which 
are not related to each other by the reciprocal lattice vectors, the triangular lattice has two 
degenerate ground states which may have some interesting consequences. For example, the 
system may be vulnerable to the Jahn- Teller distortion. 

In Fig. 4(b), we show the critical interaction U c of the Hubbard model on a triangular 
lattice as a function of total spin S z . Below U c , the BCS state ( |i~3l) is unstable against the 
paramagnetic Fermi liquid state. We find, as expected, that U c is larger than the square lat- 
tice with NNN hoping case due to the incompatibility of "pairing" and the tripartite nature 
of the lattice. Our calculated critical interaction of the 2D Hubbard model on a triangular 
lattice U c = 4.936t at S z = is a little bit smaller than the result of Krishnamurthy et. al. 
PJ20[|, their metal-insulator transition is at U C 2 = 5.271 But our result is not in contradic- 
tion with their result, because Krishnamurthy et al. have considered a more complicated 
metallic state, namely the spin density wave metallic state while we only compare our pair- 
ing state with the most simple Fermi liquid state. So our result shown in Fig. 4(b) is the 
lower bound of the critical interaction of the metal-insulator transition. 



IV. AWAY FROM HALF-FILLING 

In section II and III, we discussed the properties of Hubbard model on various lattices 
at half-filling. In HDR, this corresponds to the case where the number of doublons is equal 
to the number of holes. In this section we consider the case of Hubbard model away from 
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the half-filling. One natural choice for the wavefunction of the hole-doped system is 

n 

|^„>=E^({ft})Il7il^)>, (19) 

where 7^ = u p h K _ p + v p d p is the creation operator of one of the Bogoliubov quasiparticle 
species, |19|J21|1 \ip >= |K, S z > is the pairing wavefunction for the half-filled case (cf. Eq. 



(f[3|)). n = N h — N d , N h and N d is the number of holes and doublons, respectively. The 
wavefunction (|1^) should remain reasonably good for the multi-hole-doped system as long 
as it still has the long-range order. Especially, the wavefunction of the system with one 
extra- hole (n=l) is 

|Vi>= £4,7^0 >, (20) 
q 

which represents a series of states from the Nagaoka state to the antiferromagnetic state 
with one extra- hole as one moves from high-S^ to S z = 0. We evaluate the expectation 
value of H in the state f|20|). 



A 2 AT 

< ^\H\^ >= ]T(e fe - eK-M -N— + N d U(l + E l ({A k }). (21) 

it 

Here Ei({A k }) =< ipi\H\ipi > — < ip \H\ip > is the energy of the extra-doped hole, 

Ei(iM) = E[2A^ fc - 17(1 - 2^)v 2 k - e k v 2 k - e K W k ]\A k \ 2 - U^. (22) 

k 

From that we can find the optimum {A k } which minimizes E 1 with the constraint J2 k \Ak\ 2 — 
1 as well as the bandwidth of the doped hole which is the difference between the maximum 
and the minimum of E±. For simplicity, we consider the 2D Hubbard model on a square 
lattice with the nearest neighbor hopping only. We find that the optimum A k = 5(e^ — jJ,—rj) 
and the minimum energy of the doped hole is 

E mm = vV + -17/2 + 0*+^- f/ ^)W^ 2 + A2 - ( 23 ) 
Here rj is the optimum (e*. — /i) which minimizes the hole energy. It is given by 

rf + 77A 2 + (/i + ^- - U^)A 2 = 0, (24) 
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with the constraint 

-At- n<n <At- p.. (25) 

When S z is zero, N^/N = 1/2, /i = 0, therefore r] = 0, i.e., the hole is right on the fermi 
surface, and 



E min = A - U/2. (26) 



Here A is the solution of the gap equation (|Tj). In Fig. 5(a), we show the minimum hole 
energy E min as a function of U/t for a square lattice with the nearest-neighbor hopping only. 
One can see that in the small U region, E min decreases as U increases. It reaches a minimum 
at U/t ~ 2.7. In the large U limit, A can be written down explicitly as a Taylor series of 
t/U (cf. Eq. flT7|)). We get 

At 2 t A t A 
E mm = - ir + 0(-) = -J + 0(-). (27) 

Similarly, we find that the maximum of the hole energy at S z = is, 



E max = v/(4t) 2 + A 2 - U/2. (28) 

In the large U limit, 

Emax = 12^ + O(^) = 3J + O(^), (29) 

and the bandwidth of the doped hole, 

W = AJ + 0{^). (30) 

This is in agreement with the recent numerical calculations which find that a hole in an 
antiferromagnetic background has a bandwidth in the scale of J rather than t, in the small 
J ft limit. 122^23 

In Fig. 5(b), we plot the bandwidth 



W = E max — E min = y (At) 2 + A 2 — A (31) 
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as a function of U at S z — 0. We can see that W is proportional to t in the small U region 
where Act. It smoothly crosses over to scale as J in the large U region where A ^> t. 
Physically, this is in agreement with the string picture of a hole in an antiferromagnetic 



environment. 1^50 Due to the presence of the frustrated spin string left behind the hopping 
hole, the hole acquires a large effective mass m* which is proportional to the inverse of the 
effective hopping t*, and therefore, also to the inverse of the bandwidth W. Due to the 
presence of the interaction, t* is renormalized from t to J = At 2 /U in the large U limit. 



When S z ^ 0, the chemical potential /i is not zero anymore and is given by Eq. (|i5|). 
In Fig. 6 we show the minimum hole energy as a function of S z for U/t = 32, 16, and 8, 
respectively. We can see that E min decreases as one increases S z from zero, i.e., it costs 
less to dope a hole in the high-5 2 states than in the lower S z states. So the doping of 
holes suppresses the gap and may stabilize the high S z state. When the decrease of the 
hole energy becomes larger than the energy difference between the ferromagnetic and the 
antiferromagnetic configurations in the undoped case, the system goes through a phase 
transition from the antiferromagnetism to the Nagaoka ferromagnetism. In the large U 
limit, /i and A can be calculated analytically by Taylor expanding Eqs. ( JH] ) and flTSP, 

//= (-1/2 + N d /N)U + 0(t 2 /U), (32) 

A=[§(l-§)] 1/2 f/ + 0(t 2 /f/)- (33) 

Here N d and N is the number of doublons and lattice sites, respectively, S z = N/2 — N d . 
The S z 7^ case can again be divided into two cases due to the constraint Eq. (|25|). First, 
when At > (1/2 - N d /N)U, from Eqs. (||) and (||) we get r) = at 2 /U in the large U limit, 



here a is a number of order of 1. The minimum hole energy 

E mm = A-U/2 + 0(^). (34) 

As S z increases from zero, A decreases since the number of the doublon-hole pairs is equal 
to (N/2 — S z ). Therefore E min decrease as S z increases from zero. Second, when 4t < 
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(1/2 - N d /N)U, the optimum (t k - fi) one can get is 77 = (1/2 - N d /N)U - At due to the 
constraint Eq. fl25|), and 



E mm = vV + A 2 - Z7/2 (35) 

in the large U limit. Again, E min can be written as a Taylor series, 

27V 

E mm {S z ) = -4f(l - + 0(t 2 /U). (36) 
Finally, we have also calculated the energy of two doped holes, 

E 2 =< HH\^2 >-< ^o|#|V>o >, (37) 
using the variational wavefunction, 

|^ 2 >=E^7KVo>. (38) 

k,q 

Here A kq is determined by minimizing E 2 . The quantity E 2 — 1E\ gives information about 
the interaction between two holes. We find 

E 2 -2E 1 = — 22 v kl v k2 u k3 u k2+k3 - kl A ku k 2 +k 3 -kiAl 2k3 , (39) 

kl,k2,k3 

which goes to zero in the thermodynamic limit. This is because we did not allow the 
gap function A to fluctuate in our calculation. It is believed that the fluctuation of the 
antiferromagnetic order parameter will result in an attractive interaction between holes, 



and this attractive interaction favors a d x 2_ y 2 type of pairing. [^8| 

To summarize the results of Sec. IV, we have extended the pairing wavefunction for 
the repulsive Hubbard model to the case of less than half-filling by using the Bogoliubov 
quasi-particle operators. We use this variational wavefunction to calculate the energy of a 
doped hole. We find that in the large U limit, the doped hole has a bandwidth of order 
of J rather than t at S z — 0. |29j Here J = At 2 /U is the magnetic interaction strength. 
The hole energy decreases as one increases S z . The competition between this effect and 
the Heisenberg spin interaction which favors antiferromagnetism determines the transition 



between the antiferromagnetic state and the Nagaoka ferromagnetic state. Our variational 
wavefunction has the advantage of capable of representing a series of states in between these 
two limits in a single, relatively simple form. 

It also possible to view the system in the limit of high spin near half filling as a low 
density, two component system. Field theory methods can be applied in this limit. This 
approach lies outside the variational method of this paper, but will be described in a future 
publication. 

V. CONCLUSIONS 

In conclusion, we have studied the metal-insulator transition and magnetic ordering in 
half-filled Hubbard model using a particle-hole mapping. We show that the 2D Hubbard 
model on any lattice is an insulator at half-filling for any finite U in the (Nagaoka) high spin 
region. The physical picture of this metal-insulator transition as a function of spin is quite 
simple from the point of view of binding and unbinding of the hole-doublon pairs. At high 
spin 5* 2 , the low energy excitations are neutral, consisting of the motion of hole-doublon 
pairs. In this region, the dimensionality plays a very important role in determining whether 
the system is a metal or a Mott insulator. For d > 2, there exists a critical Ub, when and 
only when U > Ub, doublons and holes form charge neutral bound states and the system 
is a Mott insulator in the high spin region; while for d < 2, holes and doublons always 
form bound states and the system is an insulator for any finite U in the high spin region. 
As one increases the number of bound states, i.e., lowers S z , the binding costs some extra 
kinetic energy if the Fermi surface is not nested. When the cost of kinetic energy exceeds 
the binding energy, a binding- unbinding transition takes place. The excitations are charged, 
and the system is a metal. If the Fermi surface is nested, however, this binding costs no 
extra kinetic energy, and the system is always an insulator for any finite U at any spin S z . 
We have also discussed the generalization of the pairing wavefunction (|13|) to the hole doped 
case using Bogoliubov quasiparticle operators. We find that the dispersion of the doped 
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hole is significantly changed due to the presence of the interaction. In the large U limit, the 
doped hole has a bandwidth of order of J rather than t at S z — 0. 
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FIGURES 

FIG. 1. Critical interaction Ub as a function of the interlayer hopping constant t' in a 3D 
bipartite lattice. When U > Ub, a pair of hole and doublon forms charge neutral bound state and 
Hubbard model is an insulator at the high spin region. 

FIG. 2. (a). Fermi surfaces of doublons (solid lines) and holes (dashed lines) of the 2D Hubbard 
model on a square lattice with the nearest neighbor hopping only at various S z . (b). Fermi surfaces 
of doublons (solid lines) and holes (dashed lines) of the 2D Hubbard model on a square lattice with 
the nearest and the next nearest neighbor hopping at S z = 0.47V, 0.3N, and 0.2N. The parameter 
is t 2 = 0.25t. 

FIG. 3. Critical interaction U c of the 2D Hubbard model on a square lattice with NN and NNN 
hopping as a function of S z . When U > U c , holes and doublons form charge neutral bound states 
and the system is an insulator. The parameter is the same as in Fig. 2(b). 

FIG. 4. (a) Fermi surfaces of doublons (solid lines) and holes (dashed lines) of the Hubbard 
model on a 2D triangular lattice with the nearest neighbor hopping only at various S z . (b). Critical 
interaction U c of the Hubbard model on a triangular lattice with NN hopping as a function of S z . 
When U > U c , holes and doublons form charge neutral bound states and the system is an insulator. 

FIG. 5. Shows the calculated (a) minimum hole energy, and (b) the bandwidth of a single 
doped hole of the 2D Hubbard model on a square lattice with the nearest neighbor hopping as a 
function of U/t at S z = and one- hole away from half-filling. 

FIG. 6. Shows the minimum hole energy of 2D Hubbard model on a square lattice as a function 
of S z for U = 32i (solid), 16t (dashed), and 8t (dotted), respectively. 
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